Assignment

In this assignment, we are going to look at data from the Bear Brook Watershed in Maine (BBWM). The BBWM experiment has been running since the mid 1980s and is broadly focused on investigating how increased nitrogen, sulfur, and acidity affect surface waters and their surrounding watersheds.

In a previous assignment, we looked at how acid rain was changing over time and how it can control aluminum toxicity. The BBWM experiment is essentially testing what would happen if acid rain had continued unabated. BBWM is split into two watersheds, East Bear (EB) and West Bear (WB). EB is an untreated reference watershed. WB had ammonium sulfate dumped onto it on a bimonthly basis from 1989-2016. Ammonium sulfate is essentially a mixture of sulfuric acid (H2SO4) and ammonium (NH4+). This means there are a few things going on in WB:

Nitrogen availability is way up. In the 1980s the consensus was that ammonium was the most available form of nitrogen to plants.
Sulfur availability is way up. Sulfate is the most plant accessible form of sulfur. Although sulfur isn’t usually limiting, it is an important nutrient for plants and microbes.
Acidity is way up. Ammonium sulfate forms sulfuric acid in solution.

You can see the watersheds here:

'East Bear' <- st_read('data/EB/EB.shx')

## Reading layer `EB' from data source `C:\Users\gubbi\Dropbox\My PC (DESKTOP-RMMTJHU)\Documents\WR 418\WR-418-assignments\weathering_R\data\EB\EB.shx' using driver `ESRI Shapefile'
## Simple feature collection with 1 feature and 3 fields
## geometry type:  POLYGON
## dimension:      XY
## bbox:           xmin: -68.11006 ymin: 44.85966 xmax: -68.10449 ymax: 44.86458
## geographic CRS: WGS 84

'West Bear' <- st_read('data/WB/WB.shx')

## Reading layer `WB' from data source `C:\Users\gubbi\Dropbox\My PC (DESKTOP-RMMTJHU)\Documents\WR 418\WR-418-assignments\weathering_R\data\WB\WB.shx' using driver `ESRI Shapefile'
## Simple feature collection with 1 feature and 3 fields
## geometry type:  POLYGON
## dimension:      XY
## bbox:           xmin: -68.1116 ymin: 44.85935 xmax: -68.10535 ymax: 44.86438
## geographic CRS: WGS 84

mapviewOptions(basemaps = "OpenTopoMap")
mapview(`East Bear`)+
  mapview(`West Bear`, col.regions = 'red')

Q1

Part A

Given the background information above and the controls on weathering covered in lecture, write a brief (1-2 paragraphs) hypothesis about how and why you think weathering rates will differ between EB and WB.

Make sure your answer covers both what you expect to differ and why you think that will be the case.

Part B

Consider the three affects listed in the background section of this assignment (increased nitrogen, sulfur, and acid inputs). For each, give an example of where you might find this condition in the real world. Acid rain from industrial NOx and SOx has been thoroughly covered already, so pick other examples.

Q2

Data loading and manipulation.

Part A

Read in the discharge (Q) and stream chemistry data for both watersheds (see the ‘data’ folder). Combine them into a single data frame and remove the quality assurance flags from the data. Use the function glimpse() to show your final, combined data frame.

(Hint: This is the same Macrosheds data format as the Acid Rain assignment. Go look at that if you are having trouble.)

data <- read_feather('data/EB_chem_stream.feather') %>%
  rbind(read_feather('data/WB_chem_stream.feather')) %>%
  rbind(read_feather('data/WB_q.feather')) %>%
  rbind(read_feather('data/EB_q.feather')) %>%
  select(datetime, site_code, var, val)

glimpse(data)

## Rows: 227,960
## Columns: 4
## $ datetime  <dttm> 1988-03-16, 1988-03-17, 1988-03-18, 1988-03-19, 1988-03-20,~
## $ site_code <chr> "EB", "EB", "EB", "EB", "EB", "EB", "EB", "EB", "EB", "EB", ~
## $ var       <chr> "GN_Al", "GN_Al", "GN_Al", "GN_Al", "GN_Al", "GN_Al", "GN_Al~
## $ val       <dbl> 0.076, 0.082, 0.088, 0.094, 0.101, 0.107, 0.113, 0.119, 0.12~

Part B

Use unique() to show all the variables present in your combined dataset.

unique(data$var)

##  [1] "GN_Al"        "GN_ANC"       "GN_Ca"        "GN_Cl"        "GN_DIC"      
##  [6] "GN_DOC"       "GN_H"         "GN_HCO3"      "GN_K"         "GN_Mg"       
## [11] "GN_Na"        "GN_NH4_N"     "GN_NO3_N"     "GN_pH"        "GN_Si"       
## [16] "GN_SO4_S"     "GN_spCond"    "GN_TN"        "GN_TP"        "IS_discharge"

Q3

Which watershed is weathering more?

To answer this question we need to approximate weathering rates at both watersheds. We’ll be calculating annual fluxes of silicon as a proxy for chemical weathering at both sites.

Part A

Why are we using silicon fluxes and not one of the other variables you listed in the previous question?

Part B

What is the frequency of the discharge data? What about the frequency of the silicon samples? Do they match? Prove it via code. Make sure your output is able to be understood by a person.

data %>%
  filter(var == 'IS_discharge') %>%
  mutate(diff = datetime-lag(datetime)) %>%
  select(var,diff) %>%
  tail()

## # A tibble: 6 x 2
##   var          diff  
##   <chr>        <drtn>
## 1 IS_discharge 1 days
## 2 IS_discharge 1 days
## 3 IS_discharge 1 days
## 4 IS_discharge 1 days
## 5 IS_discharge 1 days
## 6 IS_discharge 1 days

data %>%
  filter(var == 'GN_Si') %>%
  mutate(diff = datetime-lag(datetime)) %>%
  select(var,diff) %>%
  tail()

## # A tibble: 6 x 2
##   var   diff  
##   <chr> <drtn>
## 1 GN_Si 1 days
## 2 GN_Si 1 days
## 3 GN_Si 1 days
## 4 GN_Si 1 days
## 5 GN_Si 1 days
## 6 GN_Si 1 days

# Yes, they match at a daily interval.

Part C

Calculate the flux of silcon at the daily timescale both EB and WB. Filter your data to only WY1998-WY2002. Use glimpse() to show your result.

The area of WB is 10.06 ha and the area of EB is 10.76 ha.

Do not remove NAs! We will need those later!

(Hint: The metadata for all Macrosheds data is available in the ‘data/macrosheds_vardata.csv’)

flux <- data %>%
  pivot_wider(id_cols = c(datetime, site_code),
              names_from = var,
              values_from = val) %>%
  select(datetime, site_code, si = GN_Si, q = IS_discharge) %>%
  filter(datetime > as.POSIXct('1998-10-01'),
         datetime < as.POSIXct('2002-09-30')) %>%
  mutate(area = ifelse(site_code == 'EB', 10.76, 10.06),
         flux = ((si*q)/area)*(86400/1)*(1/1e6)) # calc flux in mg/sec*ha and convert sec -> day and mg -> kg

Part D

Plot your new, fine-scale flux timeseries with both line and point geometry on the same plot. Be sure to include appropriate labels, a title, and a descriptive caption. Color your timeseries to match the map in the background section.

ggplot(flux, aes(x = datetime, y = flux, color = site_code))+
  geom_point()+
  geom_line()+
  scale_color_manual(values = c('blue', 'red'))+
  labs(y = 'Si Flux (kg/ha/day)', x = 'Date', color = 'Site', title = 'Daily Silicon Flux at BBWM',
       caption = 'Daily fluxes of silicon, a rock derived solute, at BBWM. \n WB has been treated with acid, EB is a reference watershed.')

## Warning: Removed 1251 rows containing missing values (geom_point).

## Warning: Removed 232 row(s) containing missing values (geom_path).

Part E

Is our data completely continuous? What will happen if we try to calculate flux with this dataset? Is there a pattern to the breaks? Why might the good folks at BBWM let those breaks happen?

(Hint: Look at timeseries of the Q and Si separately. Think the flux formula and how you would spend your research dollars.)

Part F

Let’s fill in our to make it continuous. You should be careful about doing this in your own work, but for reasons you (should) have discussed in part E we can safely do it here. There are many ways to do (with varying degrees of merit). We will be using a simple one: linear interpolation. Essentially, we are telling R to draw a straight line between observations on either side of the data gap.

To do this, we will use the na.approx() function. The ‘rule’ for linear interpolation is 2.

flux_full <- flux %>%
  arrange(datetime) %>%
  mutate(flux_int = na.approx(flux, rule = 2))

Part G

Plot your new, interpolated timeseries just like you did for part D.

p <- ggplot(flux_full, aes(x = datetime, y = flux_int, color = site_code))+
  geom_point()+
  geom_line()+
  scale_color_manual(values = c('blue', 'red'))+
  labs(y = 'Si Flux (kg/ha/day)', x = 'Date', color = 'Site', title = 'Daily Silicon Flux at BBWM',
       caption = 'Daily fluxes of silicon, a rock derived solute, at BBWM. \n WB has been treated with acid, EB is a reference watershed.')

ggplotly(p)

Part H

Zoom in on the previous gaps of your graph. Using such a naive method is… not great. But it’s fine for our purposes and it shouldn’t throw off our annual estimates by much.

We finally have what we need to calculate annual silicon flux for both watersheds. Do so and present your data as you find appropriate in units of kg per ha per year. If you color your data, make sure it sticks to the schema we started earlier.

Your figure or table should clearly answer our question: which watershed hosts the higher weathering rates?

Hints:

If you are struggling check out the work with data primer on Rstudio cloud.
The function water_year() slickly pulls the WY from a date.
Use a caption to explictly state which watershed has higher weathering rates.

flux_full %>%
  mutate(wy = water_year(datetime)) %>%
  group_by(site_code, wy) %>%
  summarize(flux = sum(flux_int)) %>%
  ggplot(., aes(x = as.character(wy), y = flux, color = site_code))+
    geom_point()+
  scale_color_manual(values = c('blue', 'red'))+
  labs(y = 'Si Flux (kg/ha/year)', x = 'Date', color = 'Site', title = 'Annual Silicon Flux at BBWM',
       caption = 'Annual fluxes of silicon, a rock derived solute, at BBWM. \n WB has been treated with acid, EB is a reference watershed. \n Note the higher fluxes at the treatment watershed.')

## `summarise()` regrouping output by 'site_code' (override with `.groups` argument)

Part I

Consider your two time series (the one with gaps and the interpolated one) and your annual estimates. Do you think your estimates are biased high, low, or are perfect? Why do you think that?

Does this affect your conclusion from the previous part? Why or why not?

Weathering

Nick Gubbins

2/21/2022

Weathering

Chemical vs mechanical

Rock derived solutes

Flux II: Return of the Flux

Does solute flux == weathering?

Importance over human timescales

Importance over geologic timescales

Assignment

Q1

Part A

Part B

Q2

Part A

Part B

Q3

Part A

Part B

Part C

Part D

Part E

Part F

Part G

Part H

Part I

Bonus